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ABSTRACT 

A finite element model for progressive in-plane damage 
and interfacial delamination was presented for cross-ply 
composite laminates. The commercial finite element 
software ABAQUS/STANDARD was used for this study. 
Each composite lamina was modeled as layers of shell 
elements based on Mindlin’s plate theory. Displacement 
continuity was established across the layers of shell 
elements. Interlaminar stresses were recovered at the 
interface with the use of multi-point constraints. Several 
material constitutive relations were developed to 
simulate in-plane damage, such as fiber breakage, 
matrix cracking and in-plane shear failure. Delamination 
was simulated by removing the displacement 
constraints. Cross-ply composite laminates were 
subjected to impact loading. The results from the finite 
element analysis compared well with the experimental 
results. 

INTRODUCTION 

The high strength-to-weight and high stiffness-to-weight 
ratios of laminated composites have made them popular 
in a wide variety of engineering applications, especially 
in high-performance structures. However, laminated 
composite structures can take maximum advantage of 
these properties only when their designs do not overlook 
the composite weaknesses due to various damage 
modes. Once the damage modes are understood, a 
better composite structure design can be achieved. 

When used in defense or automotive applications, 
laminated composites are frequently subjected to impact 
loading. They are quite susceptible to impact loading 
and their properties are severely degraded upon impact. 
Because of their heterogeneity and anisotropy within 
individual laminae and through the laminate thickness, 
they undergo complex damage process when subjected 
to impact loading. Because the damage process is 
critically important to the success of laminated 


composites, there is a great need to understand the 
details of it. 

The studies of the damage process of laminated 
composites subjected to impact loading have been 
performed both experimentally and computationally. 
Experimental studies involve preparation of test 
specimens and impact loading the specimens. Drop 
weight testing and gas gun testing are usually used for 
impact investigations. Based on the impact energy and 
absorbed energy and the inspection of the damaged 
specimens, the damage process can be studied. 

For the computational studies of the damage process, 
finite element method is frequently used. Computational 
studies provide a lot of scope for parametric studies 
leading to a better understanding of the damage 
process. But, it is very difficult to develop a 
computational model which is able to capture the 
damage process and at the same time is 

computationally efficient. This study attempts to develop 
such a finite element procedure for the study of the 
damage process of laminated composites. 

The finite element method is an effective and efficient 
tool for computational studies. Its use in fiber-reinforced 
laminated composites has been an active area of 

research for quite some time. Because of the complexity 
of the damage mechanisms that can occur in laminated 
composites, the overall mechanical behavior of 

laminated composites is highly nonlinear and is 

extremely difficult to model. In order to model the 
nonlinear response of laminated composites, continuum 
models are used. Representation of individual composite 
laminae by a homogenized continuum model permits the 
overall simulations of lamina performance and the 
damage process. 

A composite laminate can be modeled by using two- 
dimensional shell elements which use a laminate theory 
such as the classical laminate theory, Mindlin’s plate 
theory, etc. A composite laminate can also be modeled 
by using three-dimensional brick elements. Three- 
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dimensional elements are able to capture the entire 
three-dimensional stress distribution (at least 
theoretically) but they are computationally expensive. 
Thus, shell elements are the most commonly used 
elements for the study of composite laminates. A 
commercial finite element software ABAQUS was used 
for this study because user defined subroutines could be 
added to the software with ease. ABAQUS uses shell 
elements based on the first-order shear deformation 
theory in which the transverse shear strain is assumed 
to be constant through the thickness of the shell. A 
number of shell elements have been presented, based 
on higher-order shear deformation theories which could 
calculate the transverse stresses more accurately for the 
study of laminated composites [1,2]. 

The damage process in composite laminates is a three- 
dimensional process. There are intralaminar damage in 
the form of fiber damage, matrix damage and fiber- 
matrix debonding and interlaminar damage in the form of 
delamination. 

FAILURE CRITERIA 

In the analysis of composite laminates, a wide range of 
failure criteria have been used. A failure criterion may 
combine all the possible damage modes into one 
equation, e.g. Tsai-Hill failure criterion, Tsai-Wu failure 
criterion, etc. Failure criteria may also consist of a set of 
equations; each representing a distinct damage mode. 
The first two major failure criteria were proposed by 
Hashin [3, 4]; one for fiber damage and the other for 
matrix damage. The failure criteria were quadratic in 
nature since a linear criterion usually over-predicted the 
damage while a polynomial criterion was too 
complicated. The fiber damage due to tension depends 
upon the normal stress in the fiber direction and the 
shear along the fiber. The matrix damage depends upon 
the normal stress in the matrix direction and the shear 
stress along the fiber. 

Yamada and Sun [5] proposed a failure criterion of a 
lamina similar to Hashin’s fiber failure criterion but used 
in-situ strengths in the criterion. Chang and Chang [6] 
proposed phenomological failure criteria which included 
criterion for matrix cracking, fiber damage and fiber- 
matrix shear failure. The fiber damage was assumed to 
be caused by the stress in the fiber direction and the 
fiber-matrix shear failure was assumed to be due to the 
normal stress in the fiber direction and the shear stress. 
Their failure criterion for fiber-matrix shear failure was 
identical to Hashin’s failure criterion for fiber damage 
while their matrix failure criterion was also similar to 
Hashin’s. 

Chang and Lessard [7] proposed another failure criterion 
accounting for the non-linear nature of shear stress- 
strain relation in the composite materials. They also 
used in-situ strengths, similar to Yamada and Sun, in 
their criteria. Shahid and Chang [8] later on presented a 
finite element code for the analysis of composite 
materials based on this failure criterion. 


Many other failure criteria were available in the literature. 
Some of them were based on the failure of isotropic 
materials and then extended to orthotropic behavior by 
superimposing several cutoff lines that took into account 
the experimental results [9] while some used statistical 
technique to generate failure envelopes [10]. 

DAMAGE MODELS 

Once a failure criterion was met, the material underwent 
progressive damage. Different damage models were 
proposed to describe the post-failure behavior of the 
composite laminates. One of the simplest and most 
widely used methods for the post-failure modeling was 
the brittle failure model. Chang and Chang [6] proposed 
to reduce the elastic constants to zero in order to 
simulate the post-failure response of composites. This 
model was called a brittle failure model as the material 
was assumed to be perfectly brittle and there was a 
sharp drop in the strength of a damaged composite 
material. This model was used extensively for the 
analysis of composite materials [7, 8, 11, 12, 29-34]. It 
was also used for the damage in the matrix direction in 
this study. 

More recently, progressive failure models with post¬ 
failure material softening received a great deal of 
attention. In general, a macroscopic description of 
damage could be captured by a constitutive model that 
exhibited a decrease of stress with increasing strain, i.e. 
strain softening, beyond the point of failure. This strain 
softening was based on continuum damage mechanics 
(CDM) theories [13] in which the net effect of fracture 
was idealized as a degradation of the elastic moduli of 
the composite materials. All CDM models provided a 
mathematical description for the dependence of the 
elastic moduli on the damage state and the evolution of 
the damage with respect to the loading and unloading 
behavior of the composite materials. Typically, CDM 
models used a simple predefined stiffness reduction 
function, for example, a linear softening relation, an 
exponential softening relation, or more complex damage 
growth laws such as Weibull’s function. Bazant and 
coworkers presented a series of crack band models 
based on CDM [14-16]. Their models were extensively 
used for quasi-brittle materials such as fiber composites 
and concrete. 

Johnson et al. [17, 18] proposed a CDM model with 
linear softening for laminated composites. They tried to 
correlate the area under the softening zone with the 
material fracture toughness. Matzenmiller et al. [19] 
proposed a general constitutive model based on CDM 
for brittle materials. The damage function could be of 
any form. They adopted a damage function based on 
Weibull’s statistical function for illustration purpose. 
Williams et al. [20] implemented the constitutive model 
proposed by Matzenmiller in commercial finite element 
code LS-DYNA and studied impact response of 
laminated composites. They also compared this strain 
softening model with the brittle failure model proposed 
by Chang [7]. They concluded that in the post-failure 



region, the strain softening model was in a better 
agreement with the experiments when compared to the 
brittle failure model. 

Schapery and Sicking [21] used a thermodynamically 
based progressive damage formation and growth model. 
They used a polynomial evolution law for the damage in 
the matrix direction. Basu and Waas [22] extended the 
work done by Schapery and Sicking to include the fiber 
failure. They assumed that the failure in the fiber is 
brittle. Upon failure, the stress in the fiber direction 
dropped to a plateau which was maintained upon further 
loading. 

Bazant [15] showed that strain softening models were 
not objective and proposed a non-local continuum 
damage model for the analysis of brittle materials [16]. 
Kennedy and Nahan [23] implemented a non-local 
continuum damage model based on exponential 
softening in LS-DYNA for the analysis of composite 
materials. 

DELAMINATION MODELS 

The computational modeling of interlaminar failure, i.e. 
delamination, was equally challenging. A composite 
laminate is formed by stacking together individual 
laminae and these laminae can have different fiber 
orientations, resulting in non-uniform material properties 
through the laminate thickness. The mismatch of the 
material properties between the adjacent laminae could 
cause the laminae to separate from each other when 
subjected to loading. This kind of failure was called 
delamination. Liu [24] proposed a bending stiffness 
mismatch theory to predict the delamination response of 
the laminated composites. It stated that the amount of 
delamination was dependent upon the mismatch of 
bending stiffness between adjacent laminae. 
Delamination modeling posed a unique challenge in the 
composite laminate analysis as it required stress and 
displacement components in all three dimensions. 
Delamination modeling could be divided into two groups 
based on the types of failure criteria being used for the 
prediction of delamination: strength of materials based 
approach and fracture mechanics based approach. In 
both approaches the three-dimensional state of the 
composite laminate had to be known. 

Brewer and Lagace [25] proposed a quadratic criterion 
based on the interlaminar stresses for delamination 
prediction. Zhou and Sun [26] formulated the same 
interlaminar failure criterion utilizing interface stresses. 
Kim and Soni [27] proposed averaging the interlaminar 
stresses over an arbitrary critical length of one ply 
thickness. Fenske and Vizzini [28] proposed a failure 
criterion which included the in-plane stresses at the 
interface in the quadratic delamination criterion. Three- 
dimensional brick elements were used in order to 
calculate the interlaminar stresses. De Moura and 
Maques [29] used the interlaminar stresses calculated 
by a shear flexible shell element to predict the 
delamination in composite laminates subjected to impact 


loading. However, their study was limited to predicting 
the delamination shape and not to simulate 
delamination. 

Lou et al. [30] used three-dimensional brick elements 
and used Chang brittle failure criteria to study impact 
response of laminated composites. The quadratic failure 
criterion was used for delamination. Once the criterion 
was met, the element’s stress bearing capacity in the 
out-of-plane direction was reduced to zero. Other 
researchers have also used a similar finite element 
model to predict the impact response of laminated 
composites [31-34]. These models were computationally 
expensive owing to the use of three-dimensional brick 
elements. Moreover, a very refined mesh was required 
in order to calculate the inlerlaminar stresses accurately 
when using a brick element. 

Li et al. [35] proposed an interface element which joined 
the brick elements on either side of the potential 
delamination interface. Once the delamination took 
place, the stiffness of the interface element was 
changed to make the displacements discontinuous 
across the interface. Sprenger et al. [36] used a similar 
interface element to simulate delamination. 

Because of high computational cost involved in using 
brick elements, focus was shifted to model delamination 
by using computationally efficient shell elements. Multi¬ 
layered shell elements were used to model the 
composite laminates. Composite laminae on either side 
of the potential delamination surface were modeled by 
shell elements. Displacement continuity conditions were 
enforced between the shell elements. 

Zheng and Sun [37] imposed constraints between 
appropriate nodes of the shell elements to enforce 
displacement continuity. The strain energy release rate 
was calculated at the interface by the virtual crack 
closure method and the constrained conditions between 
the nodes were removed when the strain energy release 
rate was higher than the critical value. Klug and Sun [38] 
also used this model for the study of delamination. 
Wisnom and Chang [39] developed three-dimensional 
spring element to enforce the displacement continuity 
between the corresponding nodes. Their study, 
however, was based on plain strain elements. The area 
under the force-deflection curve of the spring was used 
to calculate the strain energy release rate. Once the 
delamination criterion was met, the stiffness of the 
spring was reduced to simulate displacement 
discontinuity across the interface. 

Other researchers used specially developed interface 
elements to model delamination [40-46]. Interface 
elements satisfy the displacement continuity conditions 
in a penalty sense, i.e., the stiffness of the interface 
element represented the penalty parameter. In contrast 
to continuum elements where the stress-strain relations 
were used, the relations between stresses and 
associated displacements of the interface governed the 
behavior of the interface elements. With the use of force 



and displacement components across the interface, the 
strain energy release rate across the interface could be 
calculated and the stiffness of the interface element was 
reduced to reflect the displacement discontinuity. With 
the use of interface element, the finite element modeling 
(preprocessing) was more convenient than the 
constrained conditions for the displacement continuity. 
Moreover, with the interface element, the continuity 
conditions were removed gradually, by changing the 
stiffness (penalty parameter) of the interface element 
gradually, resulting in a better representation of the 
damage process. 

IN-PLANE FAILURE 

A constitutive model is a phenomenological description 
of material behavior. In establishing a constitutive model, 
it is necessary to characterize the behavior of the 
composite material of interest by testing methods. Once 
the mechanical responses and damage modes of the 
composite material are identified, mathematical 
formulations can be established or selected to simulate 
the composite behavior. The discussions henceforth 
describe the composite responses and the incorporation 
of damage modes into a constitutive model. 

BASIC ASSUMPTIONS 

The following assumptions were made in developing a 
constitutive model for a composite lamina. 

1. A composite lamina was assumed to be 
homogeneous. Its elastic constants were averaged 
properties obtained from material tests. 
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where a is stress, e and y are strains, E is Young’s 
modulus, v is Poisson’s ratio, G is shear modulus and 
the subscripts 1 and 2 represent the directions of fiber 
and matrix, respectively. When a composite lamina 
undergoes damage, the damage may manifest itself in 
different forms. They can be divided into fiber damage 
and matrix damage. 

FIBER DAMAGE 

The stress in the fiber direction of a composite lamina is 
predominantly carried by the fibers because of their high 
stiffness in comparison with matrix material. Thus, the 
stress bearing capability in the fiber direction is not 
significantly affected by the damage in matrix. 

Fibers may fail because of tensile breakage. Broken 
fibers may debond from matrix and cause cavities. 
Fibers may also buckle or kink due to compressive 
loading. These fiber failures are catastrophic and can 
strongly affect the stress bearing capability of the 
composite lamina to which the fibers belong. The 
following failure criteria proposed by Hashin [3-4] were 
used to predict the fiber failure in this study. 


2. Each composite lamina was modeled by a two- 
dimensional orthotropic continuum. The orthotropic 
nature of the composite material was maintained 
throughout the damage process. 

3. Although the stress-strain curves for a composite 
lamina showed some degrees of nonlinearity, especially 
in shear, linear elastic relations were assumed till the 
point of failure. 


1. Fiber Failure in Tension: 0-11 -0 
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4. The damage process of a composite lamina was 
considered to be linearly brittle. This was due to the fact 
that very small permanent deformations remained in the 
composite lamina after unloading, as long as the lamina 
was not very close to complete failure. 

STIFFNESS DEGRADATION 

As mentioned above, a two-dimensional linear elastic 
orthotropic model was assumed to represent the 
behavior of a composite lamina till the point of failure. 
From Hooke’s law, the following stress-strain relation 
was given 


e 2 -(°~H)2 i j >0 failed 
- X c [< 0 elastic ^ 

X X 

where 1 and c are lamina strengths of the 
composite lamina along the fiber direction in tension and 
in compression, respectively. Fiber damage could also 
affect the integrity of matrix and hence all stiffness 
components in the constitutive model should be 
degraded when fiber damage takes place. As an 
example, Fig. 1 shows a uniaxial loading-unloading- 
reloading curve. As can be seen, it is an elastic-brittle 
response with a linear stress-strain relation till the point 
of failure. At the point of failure, the stress value drops to 
a low level and the composite lamina continues to 



maintain the same stress level on further loading. If 
unloading takes place, the material unloads with a 
secant modulus. When the stress level becomes zero, 
there is no permanent deformation in the composite 
lamina. 

A 
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where 1 is the stress to cause fiber failure. 

The degradation factor 1/1000 was chosen for the 
property degradation because if the property was 
degraded further, it gave convergence problems, while a 
degradation factor greater than 1/1000 would result in 
large post failure stress at large strains. For the 
degradation of E 1h a degradation factor 1/5 was chosen 
as it gave good match with the experimental results. 
Degradation factors of 1/3, 1/10 and 1/100 were also 
investigated. The degradation factors 1/3 and 1/10 gave 
almost similar results as the factor 1/5 but for the 
degradation factor of 1/100, the damage was too severe 
and the impactor penetrated the laminate at much lesser 
energy than in the experiments. 

MATRIX DAMAGE 


e 

Fig. 1 - Uniaxial stress-strain curve. 

When a composite lamina is damaged, the stiffness 
matrix given in Eq. (1) took the following form: 
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where the quantities in * represent the corresponding 
degraded elastic constants. They were defined as 
follows. 
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Shear stress and transverse normal stress can transfer 
to both fibers and matrix. But their effects on composite 
damage are predominantly restricted to matrix damage 
and fiber-matrix debonding. In a composite lamina, an 
initial crack in the matrix can propagate in any direction. 
However, when the crack reaches a fiber-matrix 
interface, it changes its course to propagate along the 
fiber direction without crossing into the fiber. Thus, the 
damage in matrix is aligned in the fiber direction. This 
damage of matrix depends upon the matrix strength in 
both normal direction and shearing direction, if 
debonding between fiber and matrix is not of concern. 

Being different from matrix cracking that is observed due 
to transverse tensile stress and shear stress, matrix 
crushing is observed due to transverse compressive 
stress. A number of failure criteria were proposed for 
matrix damage. In this study, failure criteria proposed by 
Hashin and modified by Chang and Lessard were used. 
It contained the following basic items. 

1. Matrix Cracking due to Shear Stress and Transverse 
Tensile Stress: °' 22 ^0 
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where { and S are the transverse tensile strength and 
shear strength of the composite lamina, respectively. 
The elastic constants were degraded to 0.1% of their 
undamaged values to simulate the brittle failure. The 
damaged stiffness matrix was defined as 
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2. Matrix crushing due to transverse compressive stress: 
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the subroutine UMAT, the material constitutive relations 
to calculate the stresses at the end of the n th increment 
and the tangent stiffness matrix for the n th increment 
were defined. 

DELAMINATION FAILURE 

Composite laminates are heterogeneous and 
anisotropic. Their properties vary from one constituent to 
another within one lamina and change from one lamina 
to another through the laminate thickness. Because of 
the mismatch in material properties, the nonuniform 
stress distribution through the laminate thickness can 
cause delamination in composite laminates. As an 
example, a [0/90] laminate has the highest mismatch of 
material properties among the [0/0] laminates, and is the 
one that is most susceptible to delamination. 
Delamination is a major form of damage as it severely 
degrades the properties of a composite laminate and 
renders the composite laminate prone to further 
damage. In fact, once delamination takes place in a 
composite laminate subjected to impact loading, it is 
very easy for the impactor to run through the composite 
laminate by cracking the matrix and pushing away the 
fibers in individual laminae. Thus, delamination is an 
important characteristic of composite damage and 
composite simulation should include delamination 
modeling. 


where c is the transverse compressive strength of the 
composite lamina. The stiffness matrix for the damaged 
composite was modified as 
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Both E22 and V 21 are defined in the same manner as 
in matrix cracking, Eqn. (7). The constitutive model 
consisting of the failure criteria and degraded stiffness 
matrices given above was incorporated into the 
commercial finite element code ABAQUS by 
programming it into a user’s material subroutine UMAT 
for the analysis of the composite material. As nonlinear 
equations are solved in increments, the solution process 
was incremental. For the n th increment, the strain 
increment for the n th increment and the stresses and 
strains of the (n-1 ) th increment were passed on to the 
subroutine by the software. The solution of nonlinear 
equations requires the tangent stiffness matrix which 
depends upon the material constitutive relations. Thus in 


DELAMINATION MODELING AND CRITERIA 

Two delamination models were commonly used, namely 
the multi-point constraint model and the interfacial layer 
model. In principle, both were based on using multi-point 
constraints to enforce the continuity of displacement 
across the laminate interfaces and to calculate 
interlaminar stresses using the constraint forces. Hence, 
both models gave similar delamination simulations. 
However, the two models differed from each other on 
how the continuity of displacements across the laminate 
interfaces was removed when the criterion for 
delamination was met. In the multi-point constraint 
model, the both rotational and translational displacement 
continuity conditions at the nodes were removed in order 
to allow displacement discontinuity. In the interfacial 
layer model, the stiffness of the interfacial element was 
reduced in order to introduce the displacements 
discontinuity. Since an additional shell element layer 
was used as the interfacial layer, the interfacial layer 
model was computationally more expensive than the 
multi-point constraint model. Thus, the multi-point 
constraint model was used in this study. 

A number of delamination models can be found in the 
literature. These criteria can be divided into two groups, 
the fracture-mechanics based models and the strength- 
of-materials based models. In the fracture-mechanics 
based models, an initial delamination front must be 
assumed along with a delamination failure criterion to 
calculate the onset of delamination. The propagation of 
delamination front, for example, can be determined by 
the strain energy release rate obtained from the virtual 



crack closure method i.e., if the strain energy release 
rate exceeds a critical level, the delamination 
propagation will take place. In the strength-of-materials 
based failure models, it is assumed that delamination 
takes place due to the presence of large interlaminar 
stresses. Failure criterion based on interlaminar stresses 
and interlaminar strengths are used to predict the 
delamination onset. In this study, strength-of-materials 
based models were used. 


enforce the continuity conditions. Accordingly the nodal 
forces in the bottom lamina will have the combined effect 
of the top and the middle laminae and thus, cannot be 
used to calculate the interlaminar stresses at the bottom 
interface. Besides, it has been noted that the bottom 
interface has the major delamination damage when 
subjected to impact loading. It was with these reasons, 
only two shell elements were used to model the three- 
lamina composite laminate. 


To begin with, a quadratic failure criterion was used to 
predict the onset of delamination. This criterion has the 
following form: 

For 0-33 -0: 
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where S is the interlaminar shear strength and ? is the 
interlaminar transverse strength in tension. It is assumed 
that delamination takes place only if the interlaminar 

normal stress 33 is in tension. This criterion was used 
by many researchers [25-34] and was also used in this 
thesis research for simulations of composite laminates 
subjected to impact loading. 

Delamination Failure Only 

A 125 mm x 125 mm composite laminate with a stacking 
sequence of [0/90/0] and thickness of 3.81 mm was 
impacted with an impactor of 12.88 kg mass at a velocity 
of 3.0 m/s. The composite laminate was modeled by two 
shell elements. The top shell element represented the 
0/90 laminae while the bottom shell element the bottom 
0 lamina. Thus, delamination was simulated to take 
place on the interface between the 0/90 lamina and the 
bottom 0 lamina. The use of a single interface in 
delamination simulations is due to the limitation of the 
commercial finite element code used in this study. For a 
laminate with two interfaces, three shell elements will 
have to be used, i.e. each shell element represents a 
lamina. However, the process of assembling the three 
laminae to become a composite laminate will render all 
the nodes but one in the laminae to become dependent 
variables. The nodes of the top lamina must be 
independent for the contact algorithm to work in 
ABAQUS Standard. Thus, the displacement continuity 
conditions are applied to the corresponding nodes of the 
middle lamina and makes the degrees of freedom of the 
nodes of the middle lamina dependent upon the nodes 
of the top lamina. Hence, the nodes of the middle lamina 
are unavailable for imposing the displacement continuity 
conditions between the middle and the bottom laminae. 
The nodes of the top lamina will have to be used to 


When the delamination failure criterion, Eq. (10), was 
used for modeling composite laminates without an in¬ 
plane failure criterion, it predicted delamination shape 
well, as shown in Fig. 2. The black elements in the 
center of the composite laminate shown in Fig. 2(a), 
collectively represent the delamination area. They were 
identified from the delaminated nodes shown in Fig. 2 
(b), to which they attached. As can be observed from 
Fig. 2, the typical peanut-shaped delamination is 
observed with the major axis aligned in the fiber 
direction of the bottom 0° lamina. 




Fig. 2 - Delamination at the bottom interface of [0/90/0] 
laminate with no in-plane failure. 

























In-Plane Failure and Delamination Failure 


When the delamination failure criterion was used in 
conjunction with an in-plane failure criterion, it was 
observed that the delamination size was greatly reduced 
as shown in Fig. 4.2. This was contrary to what was 
expected. The in-plane damage modes of the composite 
material such as matrix damage should have led to 
increased delamination. In fact, it was observed that the 
interlaminar shear stresses decreased and the 
interlaminar normal stress increased when the material 
stiffness was degraded upon matrix failure. The 
combination of the changes of interlaminar stresses 
resulted in the reduction of delamination area. Finite 
element simulations revealed that the delamination area 
based on the delamination failure criterion and an in¬ 
plane failure criterion proposed in section 3.43 was 
found to be only 2.73 cm 2 while that based on the 
delamination failure criterion only was 4.29 cm 2 . 



Fig. 3 - Delamination at the bottom interface of [0/90/0] 
laminate with in-plane failure criterion. 

Modified Delamination Failure 


It became clear that the foregoing delamination failure 
criterion must be modified. Fenske and Vizzini [28] 
proposed a delamination failure criterion with in-plane 
stress components. They claimed that delamination 
usually took place in thin matrix rich layers on laminate 
interfaces. Hence, all six stress components should be 
included in the justification of delamination onset. The 
foregoing studies seemed to support their theory. Thus, 
the delamination failure criterion was modified to include 
the in-plane stresses in this thesis research. 

The original delamination criterion proposed by Fenske 
and Vizzini takes the following form: 
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where <J I and are in-plane principal stresses, Y t 
and Y c are matrix strengths in tension and compression, 
respectively, S is interlaminar shear strength and Z t and 
Z c are interlaminar normal strengths in tension and in 
compression, respectively. 

It should be noted that both in-plane stresses and 
interlaminar stresses are included in the delamination 
failure criterion. However, Eq. (11), the distinction 
between the tensile interlaminar normal stress and the 
compressive interlaminar normal stress is not very 
strong. Although it is well known that the compressive 
interlaminar normal stress can act against delamination. 
The fact that delamination has a peanut shape 
consisting of two discontinuous lobes upon impact 
seems to support this argument. Hence, the two 
delamination lobes may not be discontinuous at the 
point of impact if the delamination failure criterion does 
not distinguish between interlaminar normal 
compression and interlaminar normal tension 
sufficiently. 

The Fenske-Vizzini delaminaiton failure criterion is 
modified to enhance the role of the compressive 
interlaminar normal stress. It is proposed that 
delamination takes place only when the interlaminar 
normal stress is in tension. Thus, the delamination 
criterion takes the following form: 


For°" 33 -0: 
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It may be noted that Eq. (12) is an extension of Eq. (10), 
with additional terms accounting for in-plane stresses. 
To calculate the in-plane stress on the laminate interface 
with the multi-point constraint model, an additional layer 
of matrix must be introduced in the finite element 
analysis. Thus, the [0/90/0] laminate which was modeled 
by two shell elements, the top element representing the 
0/90 lamina while the bottom element the bottom 0° 
lamina, is remodeled as [0/90/M/0]. The top shell 
element represented the 0/90/M lamina and the bottom 
element the bottom 0 lamina. The additional matrix layer 

























































































was very thin with a thickness of 0.01mm and did not 
have a significant effect on the global response of the 
composite laminate. With this additional matrix interfacial 
layer, it was possible to obtain the in-plane stresses on 
the laminate interface for delamination analysis. 

When the impact simulation similar to the previous 
cases was performed with the modified delamination 
failure criterion, a delamination area of 7.03 cm 2 was 
obtained. 



Fig. 4 - Delamination at the bottom interface of [0/90/0] 
laminate based on modified delamination failure 
criterion. 

SIMULATION OF [0/90/0] LAMINATES 


previous section was used to model the in-plane 
damage of the composite laminates. The delamination 
modeling presented earlier was used for modeling the 
out-of-plane damage. The impactor was modeled as an 
analytical rigid body and the shell element S4R of 
ABAQUS was used for the composite lamianae. 

Fig. 5 shows the force-deflection curves from the impact 
simulations for different impact velocities. There were 
some convergence problems in the simulations at high 
impact velocities. The simulations stopped near the 
maximum load. At High impact velocities, the 
simulations were unable to capture the penetration 
process. Figs. 6(a), 6(b) and 6(c) compare the load- 
deflection curves from those from the experiments with 
the simulations. 



- v=4.0 m/s 

-v=3.6 m/s 

- v=3.0 m/s 

v=2.7 m/s 

-v=2.4 m/s 

- v=1.9 m/s 

— v=1.7 m/s 
-v=1.2 m/s 


Fig. 5 - Force-deflection curves for [0/90/0] laminate for 
different impact velocities. 


Having proposed a constitutive model to express in¬ 
plane damage (including fiber failure and matrix 
damage) and out-of-plane damage (such as 
delamination damage), composite laminates with a 
stacking sequence of [0/90/0] were studied. Because of 
the large mismatch in material properties through the 
laminate thickness, delamination was expected to be 
severe for the cross ply laminate. Once again, 
composite laminates with dimensions of 125 mm x 125 
mm and thickness of 3.81 mm were investigated. Each 
of them was completely constrained at all four edges 
and was impacted with a spherical impactor of 12.7 mm 
diameter and 12.88 kg mass at different impact 
velocities. 



-Exp. v=3.6 m/s 

-Sim. v=3.6 m/s 


For the finite element simulations, a two-shell-element 
model was used. As the delamination model could 
handle only one interface, the delamination on the top 
interface of the 0/90 lamina was neglected, i.e. the 0/90 
laminae were modeled by one shell element. The 
second shell element modeled the bottom 0 lamina. The 
composite constitutive model of progressive fiber 
damage and brittle matrix damage presented in a 
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(a) v=3.6 m/s 



-Exp. v=2.4 m/s 

-Sim. v=2.4 m/s 


(b) 



- Exp. v=1.7 m/s 

-Sim. v= 1.7 m/s 


(c) 

Fig. 6- Comparison of force-deflection curves at three 
different impact velocities. 



(b) v=2.4 m/s 



(c) v=1.7 m/s 


The delamination shape obtained from the simulations 
for the lower interface of the [0/90/0] at three different 
impact velocities, 3.6 m/s, 2.4 m/s and 1.7m/s are shown 
in the Fig. 7(a), 7(b) and 7(c), respectively. 



Fig. 7 - Delamination shape at the bottom interface of 
[0/90/0] laminate at three different impact velocities. 

The delamination area vs. Impact energy from both 
experiments and the simulations are shown in Fig. 8. 
Clearly the simulations under predicts the delamination 
area. From Fig. 7 and 8 it can be concluded that the 
delamination model is only able to capture delamination 
damage qualitatively but is not able to match the 
experimental results quantitatively. 


















































































































































































































































































































































damage. Delamination caused an increase of about 49 
% in the energy absorption. 
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Fig. 8 - Delamination area at different impact energies. 


Fig. 9 shows the absorbed energy vs. the impact energy 
from both experimental and simulation results. The 
impact energy is based on the kinetic energy of the 
impactor while the absorbed energy is obtained from the 
integration of the force-deflection curves. 



Experiment 
Simulation 
Equal energyline 


Impact Energy, J 


Fig. 9 - Absorbed energy vs. impact energy 


To understand the effect of delamination damage on the 
response of composite laminates, two different 
simulations were performed. Both simulations were 
performed for an impact velocity of 2.4 m/s. The first 
simulation allowed only in-plane damage without 
delamination damage. The second simulation included 
both in-plane damage and delamination damage. The 
force-deflection curves from the two simulations are 
shown in Fig. 10. The energy absorbed by the 
composite laminate for the two cases were 6.41 J and 
9.91 J. The force-deflection curves and the energy 
absorbed clearly indicate the role of delamination 



Fig. 10 - Force-deflection curve for [0/90/0] laminate with 
and without delamination damage. 

When finite element simulations were performed for 
composite laminates with a stacking sequence of 
[0/30/0] and [0/45/0], there were severe convergence 
problems. Simulations with only in-plane damage and 
with only out-of-plane damage converged but 
simulations including both damage models did not 
converge. 

Fig. 11 compares the force-deflection curves of [0/90/0] 
with and without delamination, [0/45/0] with out 
delamination and [0/30/0] without delamination. All the 
simulations were performed at an impact velocity of 2.4 
m/s. It can be noted that for a hypothetical case of only 
in-plane damage and no delamination damage, the 
[0/90/0] laminate was the most rigid one followed by 
[0/45/0] and then by [0/30/0] which had maximum in¬ 
plane damage. The energy absorption was 6.41 J, 8.23 
J and 9.59 J respectively. This was due to the fact that in 
a [0/90/0], the 90-lamina arrested the cracks in the 0 
lamina due to in-plane damage to propagate. 

However, it may be noted that a [0/90/0] laminate has 
the maximum mismatch of material properties through 
the laminate thickness and thus would have the 
maximum delamination. From the earlier discussions, it 
was known that delamination degraded the laminate 
drastically. Thus, it was expected that [0/90/0] composite 
laminate would perform the worst for a real case with 
where both material in-plane damage and delamination 
damage took place. The results shown in Fig. 11 seem 
to point into this direction. Unfortunately, the simulations 
for [0/45/0] and [0/30/0] composite laminates, did not 
converge thus making the comparison with [0/90/0] 
impossible. 
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Fig. 11 - Force-deflection curves. 

CONCLUSIONS 

A finite element model was presented to study the 
damage process of composite laminates. The 
commercial finite element code ABAQUS-Standard was 
used for the computational studies. The finite element 
model used two-dimensional shell elements, thus 
making the model computationally efficient. A 
homogenized continuum model was used to 
characterize the composite laminates. Material 
constitutive relations were formulated into a user 
subroutine UMAT in ABAQUS for in-plane analysis. The 
failure criteria proposed by Hashin were used. A brittle 
damage model was used for the post- failure behavior of 
the matrix while a progressive damage model was used 
for the fiber damage. For the finite element simulation of 
the delamination damage, the multi-point constraint 
(MPC) model was used along with shell elements to 
represent composite laminate. A user subroutine MPC 
was written to impose the displacement continuity 
conditions at the interface. The nodal forces were used 
to calculate the interlaminar stresses and a stress based 
failure criterion was used for the delamination damage. 
Once the delamination was detected, the displacement 
continuity conditions were removed in the MPC model 
The delamination model was used to analyze the [0/90] 
composite beam and the results were compared with the 
elasticity analysis. The delamination shape and size 
were validated with different composite layups. A typical 
penny shaped delamination with the major axis along 
the direction of the fiber direction of the bottom lamina 
was obtained but the size predicted from the present 
model had an error of about 50 % from the experimental 
results. Two damage models, material in-plane damage 
model and the out-of-plane damage model for 
delamination, were used to investigate the response of 
[0/90/0] composite laminate subjected to impact loading 
at different impact velocities. The results from the 
simulations were compared with those from the 
experiments. The simulation results were in good 
agreement with the experiments. 
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